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1 In this short note, we discuss the basic approach to computational modeling of dynamical 

QQ ' systems. If a dynamical system contains multiple time scales, ranging from very fast to 

^v^J i slow, computational solution of the dynamical system can be very costly. By resolving the 

• , fast time scales in a short time simulation, a model for the effect of the small time scale 

variation on large time scales can be determined, making solution possible on a long 
time interval. This process of computational modeling can be completely automated. 
(N ■ Two examples are presented, including a simple model problem oscillating at a time 

scale of 10 — 9 computed over the time interval [0, 100], and a lattice consisting of large 
and small point masses. 
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1. Introduction 

We consider a dynamical system of the form 

u{t) = f(u(t),t), t£(0,T], 
t*(0) = uo, 

where u : [0, T] — > R N is the solution to be computed, uq £ M. N a given initial 
value, T > a given final time, and / : M. N x (0, T] — > M. N a given function that 
is Lipschitz-continuous in u and bounded. We consider a situation where the exact 
solution u varies on different time scales, ranging from very fast to slow. Typical 
examples include meteorological models for weather prediction, with fast time scales 
on the range of seconds and slow time scales on the range of years, protein folding 
represented by a molecular dynamics model of the form (|l.lj) . with fast time scales 
on the range of femtoseconds and slow time scales on the range of microseconds, or 
turbulent flow with a wide range of time scales. 
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To make computation feasible in a situation where computational resolution of 
the fast time scales would be prohibitive because of the small time steps required, 
the given model (|l.lj) containing the fast time scales needs to be replaced with a 
reduced model for the variation of the solution u of (|l.lj) on resolvable time scales. 
As discussed below, the key step is to correctly model the effect of the variation at 
the fast time scales on the variation on slow time scales. 

The problem of model reduction is very general and various approaches have 
been takerl^H We present below a new approach to model reduction, based on 
resolving the fast time scales in a short time simulation and determining a model 
for the effect of the small time scale variation on large time scales. This process 
of computational modeling can be completely automated and the validity of the 
reduced model can be evaluated a posteriori. 

2. A simple model problem 

We consider a simple example illustrating the basic aspects: Find u = (ui,u 2 ) : 
[0,T] -> R 2 , such that 

u x +u 1 -ul/2 = Q on(0,T], 

u 2 + ku 2 = on(0,T], (2.1) 
u(0) = (0,l) «(0) = (0,0), 

which models a moving unit point mass M\ connected through a soft spring to 
another unit point mass M 2 , with M 2 moving along a line perpendicular to the line 
of motion of M\ , see Figure [TJ The second point mass M 2 is connected to a fixed 
support through a very stiff spring with spring constant k = 10 18 and oscillates 
rapidly on a time scale of size l/v^ — 10~ 9 . The oscillation of M 2 creates a force 
~ u 2 on Mi proportional to the elongation of the spring connecting M 2 to Mi 
(neglecting terms of order u 2 ). 

The short time scale of size 10~ 9 requires time steps of size ~ 10~ 10 for full 
resolution. With T = 100, this means a total of ~ 10 12 time steps for solution of 
(|2.1[) . However, by replacing (|2.1[) with a reduced model where the fast time scale 
has been removed, it is possible to compute the (averaged) solution of (|2.ip with 
time steps of size ~ 0.1 and consequently only a total of 10 3 time steps. 

3. Taking averages to obtain the reduced model 

Having realized that point- wise resolution of the fast time scales of the exact solution 
u of (ll.l[) may sometimes be computationally very expensive or even impossible, we 
seek instead to compute a time average u of u, defined by 

i r /2 

u{t) = - u(t + s)ds, {6t/2,T-t/2, (3.1) 

T J -r/2 

where r > is the size of the average. The average u can be extended to [0,T] in 
various ways. We consider here a constant extension, i.e., we let u{t) = u(t/2) for 
* e [0,7-/2), and let u(t) = u(T - r/2) for t € (T — r/2, T]. 
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Fig. 1 . A simple mechanical system with large time scale ~ 1 and small time scale ~ 1 / yJH. 



We now seek a dynamical system satisfied by the average u by taking the average 
of (fTTTT) . We obtain 



S(f) - u(t) = /(«, •)(*) = /(u(t), t) + (/(u, •)(*) - /(«(*), *)), 



or 



u(t)=f(u(t),t)+g(u,t), (3.2) 

where the variance g(u,t) = f(u,-)(t) — f(u(t),t) accounts for the effect of small 
scales on time scales larger than r. (Note that we may extend p.2[) to (0,T) by 
defining g(u,t) = -f(u(t),t) on (0,r/2] U (T - r/2,T].) 

We now seek to model the variance g(u,t) in the form g(u,t) ~ g(u(t),t) and 
replace Q3.2j) and thus f) 1 . 1[) by 

u(t) = f(u{t),t)+g(u(t),t), te(0,T], 
5(0) = uo, 

where uo = w(0) = u(r/2). We refer to this system as the reduced model with 
subgrid model g corresponding to (jl.ip . 

To summarize, if the solution u of the full dynamical system (jl.ip is computa- 
tionally unresolvable, we aim at computing the average u of w. However, since the 
variance g in the averaged dynamical system (|3.2j) is unknown, we need to solve the 
reduced model (|3 - 3[) for u m with an approximate subgrid model g fa g. Solving 
the reduced model (|3.3|) using e.g. a Galerkin finite element method, we obtain an 
approximate solution U ~ u « u. Note that we may not expect [/ to be close to u 
point- wise in time, while we hope that U is close to u point- wise. 



4. Modeling the variance 

There are two basic approaches to the modeling of the variance g(u,t) in the form 
g(u(t), t); (i) scale-extrapolation or (ii) local resolution. In (i), a sequence of solutions 
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is computed with increasingly fine resolution, but without resolving the fastest time 
scales. A model for the effects of the fast unresolvable scales is then determined 
by extrapolation from the sequence of computed solutions^. In (ii), the approach 
followed below, the solution u is computed accurately over a short time period, 
resolving the fastest time scales. The reduced model is then obtained by computing 
the variance 

ff(u,t) = /(iv)(t) -/(«(*),*) (4-1) 

and then determining g for the remainder of the time interval such that g(u(t), t) ~ 
g(u,t). 

For the simple model problem (I2.ip . which we can write in the form (jl.ip by 
introducing the two new variables U3 = u\ and M4 = iii with 

f(u, •) = (« 3 , Ui, -Hi + ul/2, -KU 2 ), 

we note that ui ~ (for ^[kt large) while u| « 1/2. By the linearity of fx, f 2 , and 
fi, the (approximate) reduced model takes the form 

+ ux -1/4 = on(0,T], 
fi 2 + «u a = on(0,T], (4.2) 
u(0) = (0,0), fi(0) = (0,0), 

with solution u(t) = (|(1 — cosi),0). 

In general, the reduced model is constructed with subgrid model g varying on 
resolvable time scales. In the simplest case, it is enough to model g with a constant 
and repeatedly checking the validity of the model by comparing the reduced model 
(|3.3[) with the full model (|1.1[) in a short time simulation. Another possibility is to 
use a piecewise polynomial representation for the subgrid model g. 

5. Solving the reduced system 

Although the presence of small scales has been decreased in the reduced system 
(|3.3|) . the small scale variation may still be present. This is not evident in the 
reduced system (|4.2j) for the simple model problem (|2.1|) . where we made the ap- 
proximation ^2(0) = 0. In practice, however, we compute £(2(0) = - fl u%(t) dt = 
i JJ" cos(y / Kt) dt ~ 1/(i/kt) and so 112 oscillates at the fast time scale l/y/H with 
amplitude 1/( v / kt). 

To remove these oscillations, the reduced system needs to be stabilized by in- 
troducing damping of high frequencies. Following the general approactP, a least 
squares stabilization is added in the Galerkin formulation of the reduced system 
(|3.3[) in the form of a modified test function. As a result, damping is introduced for 
high frequencies without affecting low frequencies. 

Alternatively, components such as u 2 in (|4.2p may be inactivated, corresponding 
to a subgrid model of the form g2(u, •) = — •)• W e take this simple approach 
for the example problems presented below. 



March 3, 2013 22:14 WSPC/INSTRUCTION FILE 



paper 



Computational Modeling of Dynamical Systems 5 

6. Error analysis 

The validity of a proposed subgrid model may be checked a posteriori. To analyze 
the modeling error introduced by approximating the variance g with the subgrid 
model g, we introduce the dual problem 

-j>(t) = J( Ul U,t) T 0(t), t€[0,T), 

<KT) = 1>, [ ' 

where J denotes the Jacobian of the right-hand side of the dynamical system (jl.lj) 
evaluated at a mean value of the average u and the computed numerical (finite 
element) solution (/wttof the reduced system (|3.3|) . 

J(u,U,t)= [ ^-(su(t) + (l- s)U(t),t)ds, (6.2) 
Jo du 

and where ^ is initial data for the backward dual problem. 

To estimate the error e — U — u at final time, we note that e(0) = and 
4> + J{u, U, -) T </> = 0, and write 

(e(T), VO = (e(T), VO - + J(«, CT, -) T 0, e) <ft 

= r Q T (0, g - Je) eft = / T (</), 17 - it - f(U, •) + f(u, •)) A 
= Jo^, U - M ■) - g(U, ■)) dt + £{</>, g(U, ■) - g(u, •)) dt 
= R(U, •)) dt + / O T (0, g(U, •) - g(u, •)) dt. 

The first term, (c/>, R(U, ■)) dt, in this error representation corresponds to the 
discretization error U — u for the numerical solution of (13.31) . If a Galerkin finite 
element method is usecP^l, the Galerkin orthogonality expressing the orthogonality 
of the residual R(U,-) = U — f{U, ■) — g(U, ■) to a space of test functions can be 
used to subtract a test space interpolant of the dual solution <j). In the simplest 
case of the cG(l) method for a partition of the interval (0,T] into M subintervals 
Ij = (tj-i,tj], each of length kj = tj — we subtract a piecewise constant 

interpolant to obtain 

/ O T (0, R(U, •)) dt = / O T (0 - 7T0, R(U, ■)) dt < Y^ =1 ^ ma X/j \\R(U, -)\\ la J Zj \\j>\W dt 
< 5W(T)max [0)T] \\kR(U, -)\\ h , 

where the stability factor S^(T) = \\<p\\i 2 dt measures the sensitivity to dis- 
cretization errors for the given output quantity (e(T), ip). 

The second term, J (4>, g(U,-) — g(u,-))dt, in the error representation corre- 
sponds to the modeling error u — u. The sensitivity to modeling errors is measured 
by the stability factor S^(T) = \\4>\\i 2 dt. We notice in particular that if the 
stability factor (T) is of moderate size, a reduced model of the form (13.31) for 
u w u may be constructed. 

We thus obtain the error estimate 

\(e(T), VOI < (T) max \\kR(U, -)\\ h + (T) max \\g(U, •) - g(u, -)\\ h , (6.3) 
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including both discretization and modeling errors. The initial data ip for the dual 
problem (|6.1j) is chosen to reflect the desired output quantity, e.g. ip = (1, 0, . . . , 0) 
to measure the error in the first component of U. 

To estimate the modeling error, we need to estimate the quantity g—g. This esti- 
mate is obtained by repeatedly solving the full dynamical system (jl.lj) at a number 
of control points and comparing the subgrid model g with the computed variance g. 
As initial data for the full system at a control point, we take the computed solution 
U « u at the control point and add a perturbation of appropriate size, with the size 
of the perturbation chosen to reflect the initial oscillation at the fastest time scale. 

7. Numerical results 

We present numerical results for two model problems, including the simple model 
problem (|2.ip . computed with DOLFHN31 version 0.4.10. With the option automatic 
modeling set, DOLFIN automatically creates the reduced model (|3.3j) for a given 
dynamical system of the form (|1.1[) by resolving the full system in a short time 
simulation and then determining a constant subgrid model g. Components with 
constant average, such as u-^ in (|2.1j) . are automatically marked as inactive and are 
kept constant throughout the simulation. The automatic modeling implemented in 
DOLFIN is rudimentary and many improvements are possible, but it represents 
a first attempt at the automation of modeling, following the recently presented^ 
directions for the automation of computational mathematical modeling. 

7.1. The simple model problem 

The solution for the two components of the simple model problem (|2.ip is shown 
in Figure [2] for k = 10 18 and r = 10~ 7 . The value of the subgrid model g\ is 
automatically determined to 0.2495 ~ 1/4. 

7.2. A lattice with internal vibrations 

The second example is a lattice consisting of a set of p 2 large and (p — l) 2 small 
point masses connected by springs of equal stiffness k = 1, as shown in Figure [3] 
and Figured! Each large point mass is of size M = 100 and each small point mass 
is of size m — 10~ 12 , giving a large time scale of size ~ 10 and a small time scale of 
size - 10~ 6 . 

The fast oscillations of the small point masses make the initially stationary 
structure of large point masses contract. Without resolving the fast time scales and 
ignoring the subgrid model, the distance D between the lower left large point mass 
at x = (0,0) and the upper right large point mass at x = (1, 1) remains constant, 
D = V2. In Figure [51 we show the computed solution with r = 10 4 , which manages 
to correctly capture the oscillation in the diameter D of the lattice as a consequence 
of the internal vibrations at time scale 10~ 6 . 

With a constant subgrid model g as in the example, the reduced model stays 
accurate until the configuration of the lattice has changed sufficiently. When the 
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change becomes too large, the reduced model can no longer give an accurate repre- 
sentation of the full system, as shown in Figure [5] At this point, the reduced model 
needs to be reconstructed in a new short time simulation. 
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Fig. 3. Detail of the lattice. The arrows indicate the direction of vibration perpendicular to the 
springs connecting the small mass to the large masses. 




Fig. 4. Lattice consisting of p 2 large masses and (p — l) 2 small masses. 
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t t 

Fig. 5. The diameter D of the lattice as function of time on [0, 20] (left) and on [0, 100] (right) for 
m = 10~ 4 and r = 1. The solid line represents the diameter for the solution of the reduced system 
(13-31) and the dashed line represents the solution of the full system 111. 111 . 
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Fig. 6. Distance D between the lower left large mass and the upper right large mass (above) and 
the distance d between the lower left large mass and the lower left small mass (below) as function 
of time on [0, 10] and on [0, 4 • 10 — ], respectively. 



